Course: VISUAL ANALYTICS FOR POLICY AND MANAGEMENT

Prof. José Manuel Magallanes, PhD

  • Visiting Professor of Computational Policy at Evans School of Public Policy and Governance, and eScience Institute Senior Data Science Fellow, University of Washington.
  • Professor of Government and Political Methodology, Pontificia Universidad Católica del Perú.

Session 8: Spatial Data - Interactive

Let’s work with the data on contributions to Candidates and Political Committees in Washington State.

The WA portal for OpenData has this data on this website.

I have already prepare a data set, let’s get it from GitHub:

link='https://github.com/EvansDataScience/data/raw/master/contriWA.RDS'
#getting the data TABLE from the file in the cloud:
contriWA=readRDS(file=url(link))

A data frame where each row is a ZIP code, and the other columns tell us, for example, some aggregate / summary value per zip code. For example, this is the total contributed per zip code

WA_zip_contri=aggregate(data=contriWA,amount~contributor_zip, sum)

Let’s turn the numeric values into a thousand unit:

WA_zip_contri$amount=WA_zip_contri$amount/1000

Getting the Map

This is the link (please change it to yours):

mapLink="https://github.com/EvansDataScience/data/raw/master/WAzips.json"

With the help of geojsonio, we can get the map:

library(geojsonio)
## 
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
## 
##     pretty
PROJmap="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
wazipMap=topojson_read(mapLink,crs=PROJmap,stringsAsFactors = FALSE)
## Reading layer `SAEP_ZIP_Code_Tabulation_Areas' from data source `https://github.com/EvansDataScience/data/raw/master/WAzips.json' using driver `GeoJSON'
## Simple feature collection with 598 features and 102 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -124.7428 ymin: 45.54354 xmax: -116.9157 ymax: 49.0025
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Solving issues that are generally present in map files:

library(lwgeom)
## Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3
wazipMap=st_make_valid(wazipMap)
## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
library(rmapshaper)
## Registered S3 method overwritten by 'geojsonlint':
##   method         from     
##   print.location geojsonio
wazipMap=ms_simplify(wazipMap)

Merging data frame into a map:

# put map first:
layerContrib=merge(wazipMap,WA_zip_contri, 
                   by.x='ZCTA5CE10', 
                   by.y='contributor_zip',
                   all.x=F) # if no coincidence don't keep shape.

There is a new map: layerContrib.

Case: Map polygons

We will plot the the amounts contributed, which will be organised into quintiles. Let’s follow these steps:

  1. Install and load the necessary packages:
library(RColorBrewer)
library(leaflet)
  1. Decide the amount of groups:
numberOfClasses = 5
  1. Decide colors for each group of shapes (I chose a palette from here):
colorForScale='YlGnBu' # color scale
# function for COLORING quantiles in leaflet
paletteFun=colorQuantile(colorForScale,domain = NULL,
                         n = numberOfClasses)

# the base map
baseLayer = leaflet() %>%addProviderTiles("CartoDB.Positron") 
final1 = baseLayer %>%
         addPolygons(data=layerContrib,
                     weight = 0, #thickness of border
                     opacity =  1, # # the closer to 0 the more transparent
                     fillOpacity = 0.7, # color brigthness
                     fillColor = ~paletteFun(amount)) # coloring

final1

You must add a legend:

final1 %>% addLegend(data=layerContrib,
                    "bottomright",
                    pal = paletteFun, 
                    values = ~amount,
                    title = "Contributions",
                    opacity = 1)

The legend shows just percents, to get the actual intervals, you need some hard work:

final1= final1 %>% addLegend(data=layerContrib,"bottomright", pal = paletteFun, 
          values = ~amount,title = "Contributions",
          opacity = 1,
          # changes:
          labFormat = function(type="quantile", cuts, p) {
              n = length(cuts) # how many
              lower=round(cuts[-n],2) # intervals
              upper=round(cuts[-1],2)
              cuts = paste0(lower, " - ", upper) # new cuts
              }
          
     )

final1

Case: Map points:

The dataframe contriWA has columns with coordinates, which represent a point in a map. Let’s use those columns to create a spatial point data frame, while making sure it has the same coordinate system as our map:

contriWA_geo= st_as_sf(contriWA, 
                       coords = c("Lon", "Lat"),
                       crs = sp::CRS(PROJmap))

Now, plot the new map with the data from 2010:

contriWA_geo2010=contriWA_geo[contriWA_geo$election_year==2010,]

final2= leaflet(contriWA_geo2010) %>% 
        addTiles() %>%
        addCircleMarkers(clusterOptions = markerClusterOptions())
final2

Case: Several maps to higlight groups:

Imagine you need the polygons at the bottom and top deciles:

quantile(layerContrib$amount, c(.1,.9))
##        10%        90% 
##    3.44705 1906.42318

Sub maps:

#filters:
top10=quantile(layerContrib$amount, c(.9))
bot10=quantile(layerContrib$amount, c(.1))

#newMaps!
mapBot=layerContrib[layerContrib$amount<=bot10,]
mapTop=layerContrib[layerContrib$amount>=top10,]
    
legendText="Areas to watch"
shrinkLegend=0.4
title="Top and Botton Average Contribution to elections in WA (2009-2023)"

Plotting the map: And a version in leaflet:

base= leaflet() %>% addProviderTiles("CartoDB.Positron") 
final3= base %>%
        addPolygons(data=mapBot,color='blue',fillOpacity = 1,stroke = F,
                    group = "Bottom")
final3= final3%>%addPolygons(data=mapTop,color="red",fillOpacity = 1,stroke = F,
                              group = "Top")

final3

Any basic leaflet map allows interaction, but it is tricky to come back to the original situation. This is how you can do it by adding a button (check icons here:

# trick: it tell the 'center' of the state and the zoom level
textFun="function(btn, map){map.setView([47.751076, -120.740135], 7)}"


final3= final3 %>%
    
    # adding the button
    addEasyButton(
        easyButton(icon="fa-home", # a symbol
                   title="Zoom to Level 1",
                   onClick=JS(textFun)))

final3

We can use an interactive legend:

final3=final3 %>% addLayersControl(
        overlayGroups = c("Top", "Bottom"),
        options = layersControlOptions(collapsed = FALSE))
final3

Case: Arranging maps into a grid

If you need to plot more than one map:

  1. Maps will be synched:
library(leafsync)
sync(final1,final3)
  1. Maps will not be synched:
latticeView(final1,final3, ncol = 1)